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A thermodynamically-based work potential theory for modeling progressive damage 
and failure in fiber-reinforced laminates is presented. The current, multiple-internal state 
variable (ISV) formulation, enhanced Schapery theory (EST), utilizes separate IS Vs for 
modeling the effects of damage and failure. Damage is considered to be the effect of any 
structural changes in a material that manifest as pre-peak non-linearity in the stress versus 
strain response. Conversely, failure is taken to be the effect of the evolution of any mecha- 
nisms that results in post-peak strain softening. It is assumed that matrix microdamage is 
the dominant damage mechanism in continuous fiber-reinforced polymer matrix laminates, 
and its evolution is controlled with a single ISV. Three additional IS Vs are introduced to 
account for failure due to mode I transverse cracking, mode II transverse cracking, and 
mode I axial failure. Typically, failure evolution (i.e., post-peak strain softening) results 
in pathologically mesh dependent solutions within a finite element method (FEM) setting. 
Therefore, consistent character element lengths are introduced into the formulation of the 
evolution of the three failure ISVs. Using the stationarity of the total work potential with 
respect to each ISV, a set of thermodynamically consistent evolution equations for the 
ISVs is derived. The theory is implemented into commercial FEM software. Objectivity of 
total energy dissipated during the failure process, with regards to refinements in the FEM 
mesh, is demonstrated. The model is also verified against experimental results from two 
laminated, T800/3900-2 panels containing a central notch and different fiber-orientation 
stacking sequences. Global load versus displacement, global load versus local strain gage 
data, and macroscopic failure paths obtained from the models are compared to the exper- 
iments. 


I. Introduction 

Predictive capabilities of numerical tools for progressive damage and failure analysis (PDFA) of composite 
structures are predicated upon the robustness, accuracy, and objectivity of the tools. As sophisticated 
numerical tools become more feasible for widespread use in industry, significant weight and cost saving in 
the design of lightweight composite structures will be apparent. Virtual testing of materials, incorporating 
PDFA, can be used to evaluate the viability of materials or configurations prior to further scrutiny via 
physical testing. This grants analysts more flexibility during preliminary structural design stages and will 
ultimately manifest as more efficient and cost effective designs. 

Continuum damage mechanics (CDM) has emerged as a viable option for predicting the non-linear 
behavior of composite structures. The first CDM theory was developed by Refs. 1,2. Subsequently, many 
publications on this subject were produced, including numerous books. 3-7 Typically, a set of scalar damage 
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variables, or internal state variables (ISVs), introduce anisotropic damage into the composite constituent 
behavior by penalizing the components of the material stiffness tensor, and non-linear functions are used to 
control the damage evolution. Various authors have used crack density, geometry, strain energy release rate, 
and other crack features to characterize the damage evolution. 8 16 Others postulate damage evolution laws 
and characterize those laws using experiments. 17-22 CDM models must also employ failure criteria to indicate 
damage initiation. More recently, increasingly sophisticated failure criteria have been developed to better 
represent the phenomenological behavior of a damaging composite lamina 23 25 and used in conjunction with 
CDM. 

CDM techniques offer computationally efficient and readily implementable means to capturing the effects 
of damage and failure in composite materials. Unfortunately, the majority of the criteria and evolution laws 
are formulated upon phenomenological observations of the directional dependence of damage evolution, 
rather than modeling the physics of the actual damage mechanisms. Separate damage variables are used to 
degrade different components of the stiffness tensor (directions) depending on whether the damage is said 
to accrue in the matrix or fiber constituents of the composite, but the variables do not explicitly distinguish 
between the separate damage mechanisms. Furthermore, many theories involve a multitude of parameters 
that are difficult to measure and must be calibrated to correlate with experimental data. 

When implemented within the finite element method (FEM), many PDFA methodologies that utilize 
CDM breakdown when the material enters the post-peak strain softening regime locally within an element. 
Loss of positive definiteness of the tangent stiffness tensor leads to pathological mesh dependence. 26,27 
To overcome this deficiency, Ref. 28 developed the smeared crack, or crack band, model that introduces 
a characteristic element length into the formulation of the damage evolution. The original formulation 
assumed that the mode I crack band always aligns with the principle axes. 29 Ref. 30 altered the formulation 
to accommodate a fixed crack band under mixed mode conditions. An encompassing overview of smeared 
crack band models is provided by Ref. 31. 

A thermodynamically-based, work potential theory, known as Schapery theory (ST), was developed for 
modeling matrix microdamage in fiber-reinforced laminates (FRLs). 32 35 Refs. 36 and 37 extended the 
formulation to include the effects of transverse cracking by adding an additional internal state variable (ISV) 
and predicted the evolution of microdamage and transverse cracking in coupon laminates analytically. 38 
implemented this extended formulation in a numerical setting to simulate the failure of a buffer strip- 
reinforced, center-notched panel (CNP). However, due to the cumbersome nature of the evolution equations, 
the microdamage and transverse cracking evolution equations were decoupled to arrive at a more efficient 
implementation. Since no characteristic length is introduced into the formulation, the theory produces 
mesh-dependent results in a computational setting. 

The ST formulation is modified here, resulting in the enhanced Schapery theory (EST), to include the 
effects of macroscopic transverse and shear matrix cracking, as well as fiber breakage, using an approach that 
differs from Refs. 36-38. A deliberate distinction between damage and failure is made. Damage is defined as 
the effects of any structural changes resulting in a non-linear response that preserves the positive definiteness 
of the tangent stiffness tensor of the material. Conversely, failure is considered to be the consequence of 
structural changes that cause post-peak strain softening in the stress versus strain response of the material. 
Here, matrix microdamage is categorized as a damage mechanism, but macroscopic matrix cracking and fiber 
breakage are hypothesized to be failure mechanisms resulting from damage localization. The traditional ISV 
used in ST is maintained to model microdamage. Upon failure initiation, the element domain is no longer 
considered a continuum, and a smeared crack approach is used to model the embedded discontinuities. Three 
new ISVs, which incorporate the characteristic length of the finite element, dictate the evolution of the failure 
mechanisms. The EST formulation presented in Section II offers non-linear progressive damage coupled with 
mesh objective, post-peak strain softening. 

Mesh objectivity is demonstrated in Section III. In Section IV, EST is verified against experimental 
results for two center- notched panels (CNPs). Global load versus deflection data, local strain gage data, as 
well as observed failure mechanisms obtained from experiments performed at the NASA Langley Research 
Center (LaRC) and exhibited in Refs. 39 and 40 are compared to numerical results. 

II. Enhanced Schapery Theory 

The previously developed ST 34,37,38,41-43 is extended to accommodate mesh objective, post peak strain 
softening. Separate ISVs are used to govern the evolution of matrix microdamage, transverse (mode I) 
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matrix failure, shear (mode II) matrix failure, and fiber breakage (mode I). The first and second laws 
of thermodynamics are enforced, establishing thermodynamically consistent evolution laws for progressive 
matrix microdamage, as well as post-peak failure. The following sections detail the formulation of this work 
potential theory. 


II. A. Thermodynamically-based Work Potential Framework 

As a material is loaded, a measure of the work potential facilitates modeling structural changes in the 
material, such as microcracking, which affect the elastic properties of the material. Energy that is not 
dissipated is recovered when the structure is unloaded, and the magnitude of energy recovered is contingent 
upon the degraded, elastic properties at the previously attained maximum strain state. It is assumed, upon 
subsequent reloading, that the material behaves linearly, exhibiting the elastic properties observed during 
unloading, until the material reaches the preceding maximum strain state. After this state is achieved, 
structural changes resume, affecting degradation of the instantaneous elastic moduli of the material. This 
process is shown in the uniaxial stress-strain curve displayed in figure 1. The shaded area above the unloading 
line represents total dissipated potential Ws, and the triangular area underneath is the total elastic strain 
energy density We- It is assumed that the material behaves as a secant material and there is no permanent 
deformation upon unloading. This a reasonable assumption for FRLs Ref. 36; however, plastic deformation 
can also be incorporated, if necessary. 34 Extension to treat viscoelastic and viscoplastic response is outlined 
in Ref. 44. 

Both We and Ws are functions of a set of IS Vs, S m , (m = 1,2, M). These IS Vs account for any inelastic 
structural changes in the material. Differentiating Ws with respect to any ISV S m , assuming limited path- 
dependence, 34 yields the thermodynamic force, / m , available for advancing structural changes associated 
with the m th ISV. 

f = (1) 

/m dS m [ ’ 

It is shown in Refs. 33 and 34 that the total work potential is stationary with respect to each ISV. 


dW T 


= 0 


( 2 ) 


Additionally, Ref. 45 utilized the second law of thermodynamics to establish the inequality: 


fmSm > 0 


(3) 


which suggests that “healing” is not allowed for a material undergoing structural changes. Eqs. (1), (2), 
and (3) form the foundation of a thermodynamically-based work potential theory for modeling non-linear 
structural changes in a material exhibiting limited path-dependence. 


II. B. Multiple ISV Formulation of ST to Account for Multiple Damage and Failure Mecha- 

nisms 

Due to the generality of the evolution equations, Eqs. (2) and (3), the work potential theory can account 
for any number and type of structural changes that may occur in a material. This is especially useful for 
modeling progressive damage in composites because the heterogeneity of the composite, and multiaxiality of 
the local fields, enables multiple damage mechanisms to arise during a typical loading history. For instance 
in the matrix phase alone, microdamage accrues until its effects are superseded by the growth of larger 
transverse cracks. Microdamage is considered the advancement of microcracks, voids, fissures, shear bands, 
and other flaws that are present in the matrix of a composite. 36,37,42,46 The s i ze Q f these flaws is typically 
on the order of that of the fiber or smaller. Transverse cracks nucleate from pre-existing flaws within 
the matrix but grow parallel to the fibers and span the thickness of the lamina. 9,10,14, 17,47 49 Often, the 
growth of individual transverse cracks is extremely rapid; however, the effects of transverse cracking on the 
stiffness of a composite laminate can be progressive if multiple cracks form over an extended period of time 
and throughout an expansive volume. Eventually, transverse cracking is succeeded by more catastrophic 
damage mechanisms including interlaminar delamination, fiber breakage, pullout and bridging associated 
with macroscopic laminate fracture. 15,50,51 

The present EST formulation assumes that three major intralaminar mechanisms are responsible for 
all observed non-linearities in the stress-strain curve of a composite lamina: matrix microdamage, matrix 
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macroscopic cracking, and axial fiber failure. Each of these mechanisms can be accommodated by partitioning 
the total dissipated energy density, Ws, into portions associated with each mechanism. 

Matrix microdamage is the primary cause of observed non-linearity in the stress versus strain response of 
many polymer matrix composites (PMCs) (i.e. systems exhibiting negligible non-linear elasticity, plasticity 
or viscous effects) up to localization of microdamage into more severe failure mechanisms, such as trans- 
verse cracking, fiber breakage, kink band formation, or delamination. Microdamage can be considered the 
combination of matrix microcracking, micro- void growth, shear banding, and fiber-matrix debonding, figure 
2 shows a typical uniaxial response of a material exhibiting microdamage evolution, where the recoverable 
energy potential is given by W and the potential dissipated into evolving structural changes associated with 
micro damage is given by S. 

Typically, matrix microdamage continues to grow until the onset of more catastrophic failure mechanisms 
initiate. It should be noted that this work explicitly distinguishes between damage and failure in the following 
manner: 


Damage - Structural changes in a material that manifest as pre-peak non-linearity in the stress-strain 
response of the material through the degradation of the secant moduli. 

Failure - Structural changes that result from damage localization in a material and manifest as post-peak 
strain softening in the stress-strain response of the material. 

Here, three major failure mechanisms, which are distinct from the microdamage mode, are considered: 
transverse (mode I) matrix cracking, shear (mode II) matrix cracking, and axial (mode I) fiber fracture. 
These failure modes are consistent with the in-plane failure typically observed in PMC laminates. It is 
assumed that the evolution of these mechanisms yields an immediate reduction in the load-carrying capability 
of a local subvolume where the mechanism is active. Three IS Vs are used to account for mode I matrix 
cracking, mode II matrix cracking, and mode I fiber failure, respectively: S'™, Sj }, and S{. These IS Vs are 
defined completely in Section II. C, and are taken to be the potentials required to advance structural changes 
associated with these failure mechanisms. 

At any given state the total dissipated energy density Ws can be calculated as a sum of energy dissipated 
through the aforementioned damage and failure mechanisms, given by the four IS Vs. 

Ws = S + STf + S? if + S{ (4) 


According to the first law of thermodynamics, the total work potential (ignoring thermal dissipation) is given 
by the sum of the elastic strain energy density and the potentials associated with each of the damage or 
failure mechanisms. 

w T = w E + s + sy i + sfi + s{ ( 5 ) 

where We is the elastic strain energy density. Invoking the stationarity principle, Eq. (2), 

dW E 

dS 

dW E 

os? 

(6) 

dW E 

dSfi 

dW E 

dS{ 
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and the Second Law of Thermodynamics, Eq. (3), gives: 


s 

> 

0 


cm 

d IF 

> 

0 

( 7 ) 

nm 

Our 

> 

0 


s f F 

> 

0 



Eqs. ( 6 ) and (7) constitute the evolution equations for damage and failure in a material associated with 
matrix microdamage, matrix cracking, and fiber breakage in tension. 

It should be noted, that EST can also account for kink band formation under axial compression ; 41,42,52 
although, the applied loading in the examples presented in Sections III and IV are tensile, and kink banding 
does not occur. As the lamina is loaded, the fibers in the composite rotate by some angle 0, given by the 
deformation gradient in the model. To model the kink band mechanism, all calculations are then executed 
in the instantaneous fiber frame given by </>; therefore, fibers rotation induces larger shear strains, 712 . 
Increased shear strain yields more damage, leading to a reduction in the shear modulus. The increase in 
shear compliance allows for further progression of the shear strain. Under axial compression, this leads to a 
runaway instability, and a kink band will form. 


II. C. Failure Initiation 


Matrix microdamage requires no initiation criterion. For low strain levels, the microdamage ISV S remains 
small and its effects on the composite moduli are not apparent. As S evolves, with increased strains, its 
effects on the stress-strain response of the composite become more noticeable. However, it is postulated that 
the evolution of the failure mechanisms immediately yield a negative tangent stiffness; therefore, initiation 
criteria are required. Furthermore, criteria are required to mark failure initiation because the macroscopic 
cracks responsible for failure may result from localization of microdamage, or they may nucleate from pre- 
existing flaws in the material not necessarily associated with microdamage. 

EST is implemented in homogenized laminae; therefore, phenomenological criteria must be utilized that 
account for the composite microstructure. The Hashin-Rotem failure criterion incorporates separate equa- 
tions for matrix failure and fiber failure initiation. 53 The matrix failure criterion involves contributions from 
both the transverse ( 622 ) and shear ( 712 ) strains. 



€22 > 0 


622 < 0 


(8) 


where Yt is the transverse lamina failure strain in tension, Yq is the transverse failure lamina strain in 
compression, and Z is the shear failure strain. The fiber failure criterion only involves the axial strain en. 


611 

X T 


4 L ) =1 en>0 


( 9 ) 


where Xt is the maximium allowable axial strain of the lamina. A local, lamina coordinate frame is chosen 
such that, 1- is the axial direction of the fibers, 2- is the in-plane transverse direction, and 3- is the out- 
of-plane direction. When Eq. (8) is satisfied the matrix failure IS Vs S™ and Sy} are activated, and when 
Eq. (9) is satisfied fiber failure evolution Sj is permitted; otherwise, S remains the only active ISV. Upon 
satisfaction of either Eq. (8) or Eq. (9), it is assumed that the more severe failure mechanisms dominate, 
superseding the effects of matrix microdamage; therefore, 5 = 0, and additional microdamage is precluded. 
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II. D. Use of Traction-Separation Relationships to Define the Failure Potentials 

Refs. 36 and 37 used a single ISV to model the effects of transverse cracking on a composite lamina. Similar 
to microdamage, the transverse and shear moduli were related to transverse crack evolution through a set 
damage functions obtained from coupon experiments. Predictions of the non-linear response of numerous 
laminates were presented assuming a homogenous strain state in the laminates. Ref. 38 implemented 
the dual-ISV formulation of ST for predicting microdamage and transverse cracking within FEM to model 
the response of a center- notched laminate that was reinforced with buffer strips. The original formulation 
required the solution of two, coupled, bi- variate polynomials, which in an FEM framework became extremely 
computationally intensive. Thus, Ref. 38 decoupled the microdamage and transverse cracking evolution 
equations. 

In the aforementioned publications, it was assumed that the transverse cracking affected the relationships 
between stress and strain. However, the existence of a macroscopic crack invalidates the assumption of 
a continuum. Here, it is presumed that failure arises from the evolution of cohesive cracks within the 
continuum, and the IS Vs associated with failure (axial, transverse, and shear) influence the relationship 
between traction on the crack faces and the crack-tip opening displacement. The satisfaction of Eqs. (8) 
and/or (9) indicates the material behavior transitions from that of a damaging continuum to that of a 
cohesive crack, and the essential fields become traction and separation, rather than stress and strain (see 
figure 3). 

Once a cohesive crack initiates in the continuum, opening of the crack yields a reduction in traction on 
the crack faces at the crack tip. If subsequently the crack is closed, it is assumed that traction at the crack 
tip will unload linearly towards the origin of the traction versus separation law (see figure 4). The strain 
energy release rate (SERR) G J M is taken as the total energy dissipated per unit area of new surface that 
is created through crack advancement and can be calculated as the area under the traction-separation law 
(for a given traction and separation pair) minus the energy per area that can potentially be recovered by 
unloading. 

= 2 (10) 

where j indicates the material (fiber / or matrix m), M represents the corresponding mode (mode I or 
mode II), 5 j m is the crack tip opening displacement in mode M and material j, and t J M is the corresponding 
traction at the crack tip. 

Theoretically, the shape of the traction-separation laws for mode I crack growth in the fiber, and mode 

I and II crack growth in the matrix can take any shape. 54 investigated triangular, trapezoidal, beta dis- 

tribution, and sinusoidal traction-separation laws in discrete cohesive zone method (DCZM) elements and 
determined that the shape only affected convergence of the FEM solver, but not the overall results. For 
simplicity, it is assumed here that all three types of cracks obey triangular traction-separation laws, pre- 
sented in figure 4. The total area under the traction-separation curves is controlled by the corresponding 
material fracture toughness in the appropriate mode, where Gj C is the mode I fracture toughness of the 
fiber, G^c is the mode I fracture toughness of the matrix, and is the mode II fracture toughness of the 

matrix. The cohesive strengths of the materials tj c (mode I fiber strength), tf c (mode I matrix strength), 
and tf IC (mode II matrix strength) are given by the stresses in the continuum when Eqs. (8) and/or (9) are 
satisfied. Mode I, normal cracks are not allowed to grow under compression, but mode II, shear cracks can 
evolve under normal compression. Therefore, the mode I traction-separation laws for the fiber and matrix 
(figures 4a and 4b) do not accommodate negative crack tip displacements. However under negative mode 

II displacement (see figure 4c), the traction on the crack faces will increase linearly until the maximum, 
previously attained displacement magnitude is reached, after which, the crack faces will resume unloading 
according to the negative portion of the traction-separation law. The traction-separation laws exhibited in 
figure 4 do not require any initial, fictitious, pre-peak stiffness because the cracks are embedded within a 
continuum. This is an advantage over the use of DCZM elements which do require an initial stiffness because 
these interfacial elements do not actually represent physical material within the model and must attempt to 
simulate initially perfect bonding between adjacent material domains. 55,56 If set incorrectly, these fictitious 
stiffnesses can cause numerical problems. 57 

Although no mode I crack can advance under compression, it is possible for post-peak softening to occur 
under compressive loading situations. For instance a kink band could form under global axial compression, 
or the matrix could fail in local shear due to internal friction (Mohr- Coulomb) in quasi-brittle materials 
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under transverse compression . 58 Since these failure mechanisms involve local shear at a the fiber /matrix 
scales which is typically below the operating lamina/laminate scale, it appears that these mechanisms evolve 
under mode I compression. In a model containing homogenized laminae, there is no subscale shear to drive 
these compressive failure modes. However, EST could be extended further to incorporate these mechanisms 
through phenomenological accessions. The methods, developed by 42,52 and described in Section II. C, can 
be used to track the instantaneous fiber angle, and a critical fiber angle can be assigned to indicate the 
initiation of post-peak softening due to kink band formation. Similarly, a matrix compression criterion, such 
as the one developed by 23,24 could be used to signal the initiation of Mohr-Coulomb compressive failure. 
The traction-separation laws for mode I fiber compression and mode I transverse matrix compression could 
be adjusted to include the post-peak softening effects of microbuckling and Mohr-Coulomb matrix failure. 
These postulated, compressive, mode I traction-separation laws could account for energy released through 
these subscale failure modes in a homogenous model at the lamina/laminate scale. However, the examples 
presented in this chapter are tension dominated, and extension of the theory to accommodate apparent mode 
I compressive failure is left for future work. 

Using the traction-separation laws in figure 4, the SERR can be calculated with Eq. (10). 


G ? = 


t ii 


+f $f 

l IC°I 

(ii) 

+rn rra 
l IC°I 

( 12 ) 

ern :m 
l IIC°II 

(13) 

g is smeared over the entire element . 28,29 

Thus, the 


dissipation potentials in an element resulting from macroscopic cracking are related to the SERRs using the 
suitable element dimensions. 

G{ 


S{ = 


(<9+90°) 


(14) 


ST = 


am 

D II — 


qr 

7T) 

be 

G rri 

II 

iW 

be 


(15) 

(16) 


If there is a single integration point in the element, l y e ; is the length of a line running perpendicular to 
fiber direction in the element that intersects two edges of the element and the integration point, and l ^ 
is the length of a line that is parallel to the fiber direction in the element that intersects two edges of the 
element and the integration point. If there is more than one integration point in the element, the element 
can be partitioned into a number of subvolumes equal to the number of integration points, and the lengths 
^(<9+90 ) /W lengths of lines that intersect the corresponding integration point as well as two element 
edges or integration point subvolume boundaries. Incorporating a length scale into the IS Vs results in mesh 
objective, post-peak, softening. This is elaborated upon further in Section IV. 


II. E. EST Evolution Equations for a Fiber-Reinforced Lamina 

To arrive at the evolution equations for the four IS Vs, the elastic strain energy density must be defined for 
a material which may contain cohesive cracks. Therefore, the elastic strain energy We is comprised of a 
contribution from the continuum W and any possible cohesive cracks W J M . The plane stress, elastic strain 
energy density in the continuum is defined as 

w = |(£n4 + E 22 (S)e 2 22 + Gi 2 (S) 7l 2 2 ) + ^12611622 (17) 

where stress in the laminae are related to strain assuming plane stress conditions. 

&11 = Q 11^11 + Q 12^22 

&22 = Q 12^11 + Q22C22 ( 18 ) 

® Qm 712 
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where 712 is the engineering shear strain and 


C E n 

y 11 — , 

1 — ^12^21 

Q 22 = Z 

1 — Z^12^21 
Ql2 = ^12^22 

Q 66 = G \2 

V12E22 


V 21 - 


E , 


( 19 ) 


where En is the axial elastic modulus, E22 is the transverse elastic modulus, zq.2 is the Poisson’s ratio, ^21 
is the transverse Poisson’s ratio, and G12 is the elastic shear modulus. After assuming that the quantity 
^12^21 << 1 , Eqs. ( 19 ) simplify, 

Q11 = En 


Q 22 = E22 
Q 12 = V12Q22 


(20) 


Qqq — G 12 


Note that only the transverse and shear moduli (E22 and G 12) are functions of S since matrix micro- 
damage only accrues in the matrix of the laminae. The Poisson’s ratio is assumed to evolve such that the 
quantity Q12 = E 2 2^i2 remains constant; however, this restriction can be relaxed if deemed necessary. The 
degraded moduli are related to the virgin moduli (E220 and G120) and the ISV through a set of microdamage 
functions (e s (S) and g s (S)) that are obtained from three uniaxial coupon tests. 33,36,37 


E22 — E 2 2oe s (S) 


(21) 


G 12 = GMiS) 

(22) 

Degrading E22 and G12 exclusively is consistent with the intralaminar damage typically observed in PMC 
laminates. 

The elastic strain energy density of the cohesive cracks are defined as the recoverable energy per unit 
crack surface area smeared over the entire element. 

W{= ^ 

2 li e+90O) 

( 23 ) 

-l-msim 

wr = 1 1 
1 2/7 

( 24 ) 

eh 

CN 

■40 

II 

( 25 ) 

The tractions in Eqs. ( 23 )- ( 25 ) can be related to the secant stiffness’ in 

the traction-separation laws k J M . 


tj = kjSj 


tf = 

+m _ m m 
l II — K II°II 

Hence, the total elastic strain energy density in the continuum is given by 

We ( E \ ie 2 \i + ^22(5)^22 + Gi 2 ( 5 ) 7 i 2 ) + ^12^11^22 


+ 


k{(s{) f 2 kT(ST) 2 . kTiCSfi) 2 

~,(0+9O°)°I + ~,(0) + ~,(0) 


9/ (0+9O o)”I 1 JO) ^ 1 9/ (0) 

Substituting Eq. ( 29 ) into Eqs. (6) gives the ISV evolution equations. 

2 i^ 220 ^ +7 > 2Gl20 zt) = ~ 3S * 


( 26 ) 

( 27 ) 

( 28 ) 


( 29 ) 


( 30 ) 
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(31) 


2 l { e e+90O) dS{ 


dk{ gf 2 


1 dkV 


xm 2 


2/7 d5? 

1 dk 


II cm2 
°II 


= -1 


= -1 


(32) 


(33) 


2/7 dS?j 

The use of a reduced ISV S r = Si has been employed in Eq. (30). Ref. 36 has shown that the use of this 


reduced ISV yields polynomial forms of the microdamage functions in Eqs. ( 21 ) and ( 22 ). 
rule and the fact that 

r lQf ff 

aD I _ l IC 
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(34) 
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by Eqs. (11)-(16), the cohesive secant stiffnesses are 

determined. 


1 

II 


(37) 

1 

II 

r +m 

/ HdS? 

J 

(38) 

1 

II 

SS 

r +rn 

/ ,,f idS?r 
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Evaluating the integrals in Eqs. (37)- (39), while enforcing k\ 


for k J M in terms of S J M . 
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K'l ~ Z IC 
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L IC 
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results in expressions 
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L IC 


2 G 


b.m _ ,m 
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(41) 


(42) 


The thermodynamically consistent stiffnesses derived in Eqs. (40)- (42) can also be derived directly from the 
traction-separation laws using geometry. 

Finally, it is assumed that following failure initiation the strains are related to the crack tip opening 
displacements by 

^(6>+ 9 0°) eii = j(0+9O°) e C _|_ (43) 

li e h 22 = li e hg+6? (44) 

*7 712 = li e h?2 + 2 m (45) 

where ef l5 e^, and 7 ^ are the strains when Eqs. ( 8 ) and/or (9) are satisfied. Eqs. (43)- (45) imply that the 
strain in the continuum remains at the values obtained when failure initiates, and that any incremental change 
in the global strain after failure initiation is used wholly to advance the crack tip opening displacement. To 
account for changes in the continuum strain after failure initiates, it can be assumed that the stress state 
in the cracked body is homogenous and the tractions on the crack tip faces are equal to the stresses in the 
continuum. Then, the strains in Eq. (29) can be formulated in terms of the cohesive secant stiffnesses and 

the crack tip opening displacement. However, it is assumed that the evolution of strain in the continuum 
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is negligible once cohesive cracks form. Eqs. (43)- (45) can be utilized in Eqs. (40)- (42) to obtain k J M as 
functions of the global strain at an integration point. 
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Once failure initiates, the effects of failure supersede the effects of microdamage and evolution of S ceases. 
The cohesive stiffness in a cracked element is calculated using Eqs. (46)- (48) for a given strain state; then, 
Eqs. (26)- (28) and (43)- (45) are used to calculate the tractions on the crack tip faces and the crack tip 
opening displacement. It is assumed that the stress state in the integration point subvolume of the element 
is homogenous, and the tractions on the crack tip faces are equal to the stresses in the element. Lastly, the 
axial, transverse, and shear moduli of the element can be calculated: 29 
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where and G* 2 are the degraded transverse and shear moduli, due to microdamage, when Eq. (8) is 
satisfied. 

For visualization purposes in the FEM simulations, degradation parameters are defined which relate the 
current, degraded stiffnesses to their original values upon failure initiation. 
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The degradation parameter can have a minimum value of zero, which indicates that no degradation has 
occurred, or a maximum value of one, signaling that the corresponding modulus has been completely dimin- 
ished. 

The negative tangent stiffness of the stress-strain curve necessary for post-peak strain softening to occur 
imposes a restriction the maximum allowable element size, as shown by. 29 
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The analyst must be careful to ensure the dimensions of any failing elements are smaller than the conditions 
given in Eqs. (55) and (56). 

In summary, Eqs. (8) and (9) mark the transition from evolving microdamage to failure to macroscopic 
cracking. Prior to failure initiation, Eq. (30) is used to calculate the microdamage reduced ISV S r , and the 
failure ISVs S'/, S™, and S] 1 } remain zero. Eqs. (21) and (22) are used to calculate the degraded transverse 
and shear moduli. Subsequent to matrix failure initiation, microdamage growth is precluded, and S r remains 
at S*, the value of S r when Eq. (8) was satisfied. The degeneration of the transverse and shear moduli, 
resulting from matrix transverse and shear cracking, is calculated using Eqs. (50) and (51). Finally if Eq. 
(9) is satisfied, the axial modulus is calculated using Eqs. (49) as fiber breakage evolves in the element. 
Once the material moduli have been calculated using the appropriate evolution equations, the stresses can 
be updated accordingly using Eqs. (20). 


III. Mesh Objectivity 

The theory outlined in Section II eliminates the mesh dependency that arises from ill-posedness when the 
elements enter the post-peak softening regime by introducing a characteristic length into the formulation. 
The total SERR dissipated during the evolution of the discontinuity is equal to the prescribed fracture 
toughness and is independent of the element size. This approach, commonly referred to as the smeared 
crack approach, or crack band model, has been used to alleviate mesh dependency in FEM since it was first 
developed by 28 for post-peak strain softening in concrete. 

To exhibit the mesh objectivity of EST, a simple example is presented in this section. One quarter of a 
200 mm x 100 mm panel containing a central hole is modeled with finite elements using the Abaqus, version 
6.10-1 finite element software. 59 The panel contains a hole with a radius of = 5 mm in the center. The 
left edge of the panel is constrained in the x-direction to simulate symmetry. Similarly, the bottom edge is 
constrained in the ^/-direction. A uniform displacement is applied to all the nodes on the top edge of the 
panel in ^/-direction. Details on the panel geometry and boundary conditions are displayed in figure 5. The 
panel is composed of a generic, [90°], composite lamina with the fiber angle measured with respect to the 
y- axis; thus, the applied displacement is perpendicular to the fiber direction in the panel. EST is used to 
model damage and failure in the panel. 

Four different meshes are used to evaluate the effect of mesh size on the response of the panel. All four 
meshes consist of two-dimensional (2-D), plane stress, quadrilateral, S4R shell elements. The elements are 
linear, reduced integration elements and contain four nodes and one integration point each. The density 
of the four meshes, shown in figure 6, increases within a region near the central hole. Average element 
sizes equal 0.5^, 0.2 r^, O.lr^, and 0.04^ are used in the four different meshes. Coarser elements are used 
away from the hole to improve computation time. The four meshes are subjected to the same boundary 
conditions and loading, and are composed of the same material properties. The same elastic, damage and 
failure parameters are also used in all four simulations. 

The resultant, applied tensile stress (given by two times the sum of the reaction forces at the nodes on 
the top edge divided by the cross-sectional area) normalized by the transverse, mode I, critical strain times 
the transverse Young’s modulus a is plotted versus the applied displacement normalized by the radius of the 
hole A for the four different meshes in figure 7. It can be seen that the mesh density has a minimal effect on 
the load-deflection results. The small discrepancy in the results between the four meshes can be attributed to 
the increased accuracy in the fields as the mesh is refined. Moreover, the total energy dissipated is preserved 
from mesh to mesh. 

Figure 8 displays contours of the normalized, reduced microdamage ISV S r immediately before failure 
initiation in the four different meshes. S r is normalized by the maximum S r obtained in all four simulations 
which is 27% of the maximum allowable S r required to bring the moduli to zero. The four meshes exhibit 
similar microdamage contours, but as the mesh is refined, the microdamage is contained to the vicinity of 
the hole. Additionally, the global stress at which failure initiates reduces slightly as the mesh is refined; this 
supports the previous statement that the discrepancies in the stress-displacement responses were a facet of 
increasing field accuracy with mesh refinement. 

The transverse matrix failure degradation parameter DJ 1 is plotted for all four meshes in figure 9. figures 
9a-9d show the failure pattern at the ultimate, global load. The coarsest mesh shows that a crack band 
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has grown at the intersection of the hole and the bottom symmetric boundary and is propagating towards 
the free edge, while moving away from the bottom boundary when the specimen ultimate load is achieved. 

In the finer meshes, the crack band path is different. In figures 9b-9c the crack band initiates at the hole 
slightly above the bottom boundary and propagates towards the free edge while approaching the bottom 
boundary. Since crack band initiation is governed by the quadratic H-R failure criterion (Eq. (8)), both the 
transverse and shear strains are contributing to the initiation leading to crack band initiation at the hole 
boundary away from the symmetric, horizontal edge where the transverse strain is not the maximum. Any 
failure criterion, in theory, could be used to govern crack band initiation, and if a maximum strain criterion 
was used, the crack band would always initiate at the hole boundary at the symmetric edge of the model, 
which is the location of the maximum transverse strain. The finest mesh, figure 9d, exhibits multiple crack 
bands near the hole, but only one of the crack bands grows significantly. While in all simulations the first 
crack band initiates at the hole, figures 9c and 9d show some crack bands beginning to initiate at the free 
edge and propagate inwards by the time the ultimate load is attained. Figures 9e-9h display the transverse 
crack band pattern once the specimen has lost all load-carrying capability. The solution for the simulation 
with the finest mesh 0.04?^ diverged before all the load carrying capability was lost; so, figure 9h presents the 
crack band path at the final converged state, which is still far below the ultimate load state. In figures 9e-9h 
the same crack band patterns that developed in figures 9a-9d are evident, except those crack bands have 
saturated to maximum degradation: D™ = 1. Additionally, the crack bands advancing from the free edge 
observed in figures 9c and 9d have progressed further. However, this is well after the specimens reached their 
ultimate load; therefore, it is assumed that those crack bands did not influence the load carrying capability 
of the panels. The discrepancy in crack band path observed for the different meshes indicate that the crack 
band path is dictated by the accuracy of the fields surrounding the leading boundary of the crack band path. 

IV. Example - Center Notched Panels Subjected to Uniaxial Tension 

IV. A. Experimental Details 

Two center- notched panel (CNP) configurations were tested at the NASA Langley Research Center (LaRC). 39,40 
The geometrical details of the panel and testing boundary conditions are presented in figure 10. The panels 
were 3” wide and 11.5” long. Two 3” x 2.75” tabs were placed on both ends of the specimens, leaving a gage 
section of 3” x 6” which is displayed in figure 10. A central notch was machined in each panel that was 0.75” 
wide and had a notch tip radius of 0.09375”. The end tabs were clamped and a vertical, tensile displacement 
(in the //-direction) was applied to the top tab using a servo-hydraulic testing machine. The bottom tab 
was fixed preventing any //-displacement of the bottom boundary of the gage section. The gripped tabs also 
prevented any displacement in the x-direction at the top and bottom boundaries of the gage section. 

The panels were comprised of laminated T800/3900-2 carbon fiber /toughened epoxy composites. Three 
different lay-up configurations were tested; however one of the configurations exhibited significant delamina- 
tion. Since the focus of this work is modeling in-plane damage and failure mechanisms, this configuration is 
not considered here. The two remaining configurations are presented in table 1. The first lay-up, Laminate- 1, 
consists of 12 0° plies, and the second, Laminate-2, is a symmetric, multi-angle lay-up with 40% 145° I, 40% 
0°, and 20% 90° layers. 

Several strain gages where affixed to the test panel, labeled Sg-1 through Sg-4 in figure 10. Sg-1 was 
placed in the center of the panel, 1.5” above the notch. Sg-2 was placed 1.5” above the notch tip. Sg-3 was 
attached in front of the notch, 0.5” from the free edge, and Sg-4 placed at the notch tip. Global load versus 
displacement data, and local strain gage data was reported in Refs., 39,40 along with a post-test C-Scan of 
Laminate- 1 and photograph of Laminate- 1 and photograph of Laminate-2. 

IV.B. Finite Element Model Details 

The linear elastic properties of T800/3900-2 used in the FEM models are presented in table 2, and were 
taken from Ref. 39. The shear microdamage function g s utilized in Eq. (22) was obtained from [45° /-45°]ss 
angle-ply T800/3900-2 coupon tests. The transverse, tensile and compressive microdamage functions were 
inferred by scaling the coefficients of the microdamage curves presented by 36 for AS4/3502 by the ratio of 
the virgin transverse modulus of T800/3900-2 to that of AS4/3502, as the stress-strain curves of the coupon 
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laminates necessary to characterize e s were not available. The polynomial forms of e s and g s are 

e s (S r ) = e S Q + e s i S r + e S 2 + e s sSf + e S 4 $jf + (57) 

9s{S r ) = #s0 + + 9s2$r + 9s3^f + ^4^ + Qs^Sr (58) 

The coefficients of the microdamage curves are presented in table 3, and the curves are plotted in figure 12. 

To increase computational efficiency, the first derivatives of the higher order microdamage polynomials 
are approximated using linear spline interpolants. Thus, the microdamage evolution equation, Eq. (30), is 
always second order in S r yielding a very efficient analytical solution. Since the value of S r from the previous 
increment is used to estimate which spline regime should be used to solve for S r in the current increment, 
the solution is checked to ensure that it falls within the applicable range of the spline that was used. If it 
falls outside of the range of S r that are valid for the splines, the solution is calculated again using splines 
that accord to the solution of the previous iteration. This procedure continues until the solution of Eq. (30) 
falls within the relevant range of S r for the splines used in Eq. (30). A maximum number of iterations 
can be set, after which Eq. (30) is solved using the full polynomial forms of the damage functions, given in 
Eqs. (57) and (58), by finding the eigenvalues of the companion matrix of the polynomial coefficients of the 
evolution equation. 

The axial mode I, transverse mode I, and shear mode II critical cohesive strains, and fracture toughness’ 
are given in table 4. The matrix mode I and mode II cohesive critical strains ( Y T , Yc, and Z) and the fracture 
toughnesses (G^c and G^} c ) were calibrated using data from Laminate- 1. These values were adjusted until 
the global load versus local strain at strain gage Sg-1 (see figure 10) obtained from the model provided the 
best match to the experimental data. Both the simulation results and experimental data are presented in 
figure 14a. Laminate- 1 did not exhibit any axial failure; so, the fiber mode I parameters (Xt and Gj C ) were 
calibrated such that the ultimate load from the simulation of Laminate-2 corresponded with the ultimate 
load reported by 39 for Laminate-2. Subsequent work will outline a procedure for measuring these values 
experimentally. 

Displacement was applied to both laminates using the *DYNAMIC keyword in Abaqus with the parame- 
ter APPLICATION = QUASI-STATIC. This implicit dynamic solver is recommended for quasi-static problems 
exhibiting a high-degree of nonlinearity. This procedure uses numerical damping to stabilize the problem. 
The numerical damping does not significantly affect the simulation results because the velocities in these 
simulations are low. For Laminate-1 a total displacement of 0.0236”, and for Laminate-2 a displacement of 
0.0472”, is applied over 1000 seconds. The panels were assigned a representative density of 0.057 lb/in. 3 . 
This implicit, dynamic technique has advantages over traditional static, implicit solvers which have difficulty 
converging when the material exhibits post-peak softening, 60,61 and is not limited by a minimum stable time 
step required with explicit solvers. 62 

IV. C. Results - Laminate-1 

Global load P versus displacement A of a 4” section of Laminate- 1 is compared to results from the EST 
simulation in figure 13. Very good agreement between the model and the experimental results is achieved. 
The response of the specimen appears to be linear until near 8,000 lbf., where the specimen begins deforming 
non-linearly. The EST simulation captures the initiation and progression of the global nonlinearity accurately. 
This panel was not loaded until catastrophic failure; hence, the data presented in figure 13 represents load 
versus displacement data prior to the ultimate load of the specimen. 

Local strain gage data (global load P versus local ^/-direction strain e yy ) from Laminate- 1 is plotted with 
the results from the EST FEM model in figure 14; please refer to figure 10 for locations of strain gages. 
Strain relaxation is observed in the gage farthest away from the notch: Sg-1, shown in figure 14a. The 
mode I and mode II matrix failure parameters in EST were calibrated such that the model demonstrates 
the same transition into strain relaxation at this location and at a similar global applied load. This load, 
taken as the splitting load, is 8,250 lbf. in experiment and 8,210 lbf. in the model (summarized in table 
5). The transition to strain relaxation is more abrupt in the experiment as evidenced by the sharp knee 
in the load-strain curve, whereas, the transition in the model is more gradual. Prior to strain relaxation 
at this point, the experiment displayed slight stiffening not observed in the model. Additionally, the model 
response is much smoother than that of the experiment in the strain relaxation regime. Even though the 
global loading is quasi-static, local events, such as cracking, may be dynamic; therefore, the discrepancy in 
the strain relaxation portion of the load-strain curves could be a result of dynamic matrix crack growth and 
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arrest in the test specimen. Local crack dynamics were not taken into account in the model. Additionally, 
the jaggedness of the experimental data may be a facet of the stochastics related to the local microstructure 
of the composite that are not included in the model. The data from the experiment and simulation for Sg-2, 
which is located 1.5” directly above the notch, are presented in figure 14b. The model predicts less strain at 
Sg-2, for a given load, than the experiment, but the non-linear trends are very similar. This gage lies directly 
in front of the splitting crack path, shown in figure 10, and it is not realistic to expect perfect agreement in 
areas experiencing high levels of damage and failure because of idealizations used to model the evolution of 
cracks in the simulations. Figure 14c displays data for Sg-3, located in front of the notch near the free edge. 
Very good agreement between the experimental and simulation results are exhibited. The model accurately 
captures the non-linear evolution of strain, away from the highly damaged regions, as a function of applied 
load. Finally, results for Sg-4 (located directly at the notch tip) are given in figure 14d and includes the 
experiment and simulation display of axial strain relaxation. As with Sg-2, Sg-4 shows less strain for a 
given applied load. However, the load at which the strain at Sg-4 relaxes in both the experiment and model 
correlate well, in accordance with the splitting load. Again, the relaxation response of the experiment is 
discontinuous, but the model exhibits continuous behavior. 

A C-Scan of the failed Laminate- 1 specimen is displayed in figure 15. Four splitting cracks can be 
observed propagating outward from the notch tip, parallel to the loading direction, towards the gripped 
edges. Contour plots of the normalized microdamage obtained from the simulation are presented in figure 16 
at the splitting load 8,200 lbf. and at 16,400 lbf. In these plots, S r is normalized by the maximum achievable 
value, S™ ax = 7.57 psis, obtained from figure 12. In Laminate-1, S r reached a maximum value equal to 
O.UIS™ 0,0 ^ . At the splitting load, the regions of maximum damage are localized to small regions, along the 
same crack path observed in figure 15, embedded in a more widespread domain and exhibiting less severe 
microdamage. A similar microdamage contour is observed at P — 16,400 lbf., except the localized damage 
region has nearly proceeded to the fixed boundaries of the panel. 

Figure 17 shows the shear failure degradation factor D j} at the splitting load and 16,400 lbf. The shear 
failure localizes into crack bands that are a single element wide and progress equivalently to the cracks 
observed in the experiment. At the splitting load (figure 17a), the crack bands have progressed less than half 
of the way between the notch and the panel boundary on either side of the notch. In figure 17b, the crack 
bands have nearly reached the fixed grip boundaries. Additionally, figure 17b displays some irregularity in 
the crack path. In these regions, the mesh is not uniform because it needed to accommodate larger elements 
used to represent the strain gages (see figure 11). No axial failure was observed in this simulation. The shear 
crack path was in a state of transverse compressions, and according to the traction-separation relations used 
in figure 4b, transverse failure does not progress under compressive conditions. Thus, contours of and 
D™ are not shown. 

IV. D. Results - Laminate-2 

Numerical results for applied load versus displacement of a 4” section of Laminate-2 are presented in figure 
18. The experimental ultimate load 15,300 lbf. correlates well (axial failure parameters were calibrated to 
obtain an ultimate load that most closely matched the experimental data) with the ultimate load obtained 
from the model, also 15,300 lbf., and is summarized with the splitting load from the Laminate- 1 analysis in 
table 5. The global response up to failure is nearly linear and failure occurs suddenly and catastrophically. 

Figure 19 compares the applied load versus strain gage results from the model to the data from the 
experiment. Sg-1 and Sg-2 exhibited similar behavior; the strain increases until the ultimate load is obtained, 
after which the strain relaxes abruptly. The experimental data and numerical results both display this 
behavior. The model exhibits slightly more strain, for a given load, prior to ultimate failure. At Sg-3, the 
model predicts strain localization after the ultimate load is achieved. The gage data shows a slight reduction 
in strain as the load drops; however, the gage was placed directly in the crack path and may have been 
damaged when the panel failed. The model results and experimental data for Sg-4 exhibit similar trends, but 
the strain gage shows a large degree of nonlinearity at the notch tip. 39 attributed this observed nonlinearity 
to local interlaminar stresses near the notch free edge which caused some local delaminations. Since the 
focus of this work was modeling in-plane damage mechanisms, these effects are not captured; however, the 
model could be easily extended to incorporate delamination by placing DCZM elements between continuum 
shell layers. 40 

In the experiment, the gages measure the strain over a continuous area associated with the size of the 
gage, but in the model, the strain is taken at the integration point of an element; thus, these measures should 
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not be expected to correspond exactly. In areas where there are large gradients present, such as near a notch 
tip (Sg-4) or near cracks (Laminate-1, Sg-2 or Laminate-2, Sg-3), it becomes even more difficult to relate 
the strain gage data to numerical strains from a discretized continuum. This may contribute to some of the 
discrepancies between the local gage data and the model results in figures 14 and 19. 

A photograph taken of the failed, Laminate-2 specimen is presented in figure 20. The photograph shows 
that two macroscopic cracks initially propagate from the notch tip towards the free edges, perpendicular to 
the applied load, in a self-similar fashion. Eventually, the cracks turn and proceed towards the free edge at 
an angle. 39 claim, supported by visual image correlation displacement data, that there was some eccentricity 
in the specimen alignment, which resulted in deviation from self-similar crack growth. 

Normalized microdamage contours just prior to the ultimate load are presented for the outermost 45°, 
0°, -45°, and 90° plies in figure 21. Similar microdamage patterns are evident in the 45° and -45° layers. 
Microdamage propagates outward, toward the free edge, from the notch tip in petal-like patterns. The 
microdamage in these layers is highly distributed throughout the plies. The 0° ply displays a more contained 
microdamage pattern along the lines of the microdamage contours associated with axial splitting as shown 
in figure 16. A moderate level of microdamage is also displayed in the 90° layer, but a low degree of 
microdamage is distributed throughout most of the layer. 

Figure 22 shows the axial failure degradation parameter Dj at the ultimate load for the four unique layers. 
A small amount of axial failure in the 45°, 0°, and -45° layers can be observed at the notch tips. It appears 
that more failure occurs at one notch tip than the other. This can be attributed to numerical imperfections 
resulting from dissimilar meshes at the opposite notch tips, that is the mesh is not symmetric about the 
y- axis. No axial failure is observed in the 90° layer. Contours of the transverse, mode I, failure degradation 
parameter D J 1 at the ultimate load are plotted in figure 23. The failure patterns are similar in the 45° and 
-45° plies in figures 23a and 23c and are comparable to the microdamage contours in figures 21a and 21c, 
except the failure is restricted to regions on either side of the notch. Furthermore, small, highly degraded 
domains can be observed propagating from the notch tip at an angle corresponding to the fiber direction in 
the ply. The 90° layer exhibits some moderate degradation in a localized region around the notch tips, and 
the 0° layer does not exhibit much D™. Contours of the shear, mode II, failure degradation parameter DJ} 
are presented at the ultimate load in figure 24. Very similar failure paths can be seen in the 45° and -45° 
layers and the patterns are nearly symmetric across both centerlines of the panel. This is expected because 
as figure 4c indicates, the sign of the local shear strain does not affect the failure degradation. Dy\ in the 
0° and 90° is limited to very small regions surrounding the notch tips. 

Contours representing the microdamage in the four unique layers are presented in figure 25 after the 
panel has completely failed and lost all of its load carrying capability. Although further matrix microdamage 
evolution is prohibited in elements that have failed (transverse/shear or axial), in the other elements that 
have not failed, matrix microdamage evolution continues. Nearly the entire 45° and -45° layers reach a 
microdamage level of 0.1 8S™ ax . The 0° and 90° plies exhibit similar microdamage patterns; however, low 
levels of microdamage are more widespread in the 90° ply. figure 26 shows the fiber failure path once the 
specimen has completely failed. All of the layers, except the 90° layer, show self similar cracks propagating 
from the notch tips towards the free edges of the panel. The angled crack path shown in figure 20 was not 
reproduced because the eccentric loading (suspected in the test) was not introduced into the simulation; 
therefore, the crack growth remained self-similar. A high degree of transverse matrix failure can be seen 
in the axial crack path in the 45°, -45°, and 90° plies in figure 27. In the 0° layer, some transverse failure 
is observed surrounding the fiber failure, as well as away from the axial failure path, which resulted from 
a stress wave reflecting off the free edges when the axial crack band reaches the boundary. Finally, D j} is 
presented after the specimen has failed in figure 28. Similar failure to figures 24a and 24c in the 45° and -45° 
is exhibited, but a highly degraded region has localized in the axial crack path, figures 28b and 28d show 
fairly extensive regions containing a high degree of shear matrix failure surrounding the axial failure path. 

V. Conclusions 

A thermodynamically-based, work potential theory for damage and failure in composite materials, en- 
hanced Schapery theory (EST), was developed. A marked distinction between damage and failure was 
introduced. Damage was considered to be the evolution of mechanisms that cause structural changes in the 
material such that the non-linear tangent stiffness tensor remains positive definite. Failure was taken to 
be the effect of structural changes in the material that result in loss of positive definiteness of the tangent 
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stiffness matrix and post-peak strain softening. Separate internal state variables (IS Vs) were used to account 
for damage and three in-plane failure mechanisms. 

In EST, matrix microdamage, which includes matrix microcracking, shear banding, and microvoid growth, 
is responsible for all damage in a composite lamina and was accounted for with a single ISV, along the lines 
of the original Schapery theory (ST) formulation. The relationship between the transverse and shear moduli 
of the lamina were related to the ISV through a pair of experimentally-obtainable microdamage functions. 

Three major, in-plane failure mechanisms applicable to continuous fiber-reinforced, laminated, polymer 
matrix composites were identified: mode I matrix cracks, mode II matrix cracks, and fiber breakage. A failure 
initiation criterion was used to mark the transition from a damaging continuum to a damaged continuum 
with an embedded discontinuity. After failure initiation, microdamage evolution ceases and separate IS Vs 
are introduced to incorporate the effects of the three major failure mechanisms. Evolution of the failure IS Vs 
is based upon traction-separation laws (which are a functions of the appropriate fracture toughnesses) and a 
characteristic element length. Typically, the existence of a non-positive definite stiffness tensor would result 
in pathologically mesh dependent solutions; however, in EST, mesh objectivity is ensured by incorporating 
a characteristic length scale into the failure evolution. 

In Section III of this paper, mesh objectivity is demonstrated. A unidirectional composite plate with a 
central hole obeying EST is loaded in transverse tension and the response is calculated using four different, 
increasingly refined, meshes. The global stress versus strain response remained unaffected by the change 
in mesh, save for the effects from increased accuracy of local fields in the vicinity of the hole with denser 
meshing. 

Two center- notched panels, with different lay-ups, composed of T800/3900-2 were tested under tensile 
loading at NASA LaRC. Global load versus displacement and global load versus local strain gage strain 
data were compared to results obtained from FEM models utilizing EST in Section IV. Quantitatively, very 
good correlation was achieved for both laminates. Furthermore, damage and failure paths predicted by the 
models matched well with the experimental results. 
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ID 

Stacking Sequence 

Thickness (in.) 

Laminate- 1 

[0°]l2 

0.078 

Laminate- 2 

[45707-45 o /0790 o ] S 

0.065 


Table 1: T800/3900-2 lay-up configurations used in CNP tests at NASA LaRC. 


Property 

T800/3900-2 

En (Msi) 

23.2 

E 22 (Msi) 

1.3 

G \2 (Msi) 

0.9 

M2 

0.28 


Table 2: Linear elastic properties for T800/3900-2 used in FEM models. 


e s (S r ) Coefficients 

Tension 

Compression 

g s {S r ) Coefficients 



1.0000 

1.0000 

9s0 

1.0000 

e s i 

-6.0373E-2 

8.4887E-4 

Qsi 

-3.0567E-2 

e S 2 

2.5937E-2 

2.8002E-2 

9s2 

-1.2135E-1 

&s3 

-1.5789E-2 

-6.2122E-3 

9s3 

3.7438E-2 

&s4 

2.2571E-3 

N/A 

9s4 

-4.5405E-4 

5 

-1.0440E-4 

N/A 

9s5 

1.9532E-4 


Table 3: Microdamage function coefficients for T800/3900-2 used in FEM models. 
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Property 

T800/3900-2 

X T 

0.021 

Y t 

0.0092 

Y c 

0.0115 

z 

0.0075 

nf 

U IC 

1026 lb £T 

U IC 

2.39 

rrm 
U IIC 

6.78 


Table 4: Failure parameters for T800/3900-2. 



Type 

Experimental 

Numerical 

Laminate- 1 

Splitting 

8,250 lbf. 

8,210 lbf. 

Laminate- 2 

Ultimate 

15,300 lbf. 

15,300 lbf. 


Table 5: Critical experimental and simulation loads for Laminate- 1 and Laminate-2. 
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o 



Figure 1: Typical stress-strain curve, containing pre-peak nonlinearity and post-peak strain softening, 

showing the total elastic (We) and total dissipated (Ws) potentials. 



Figure 2: Typical stress-strain curve with a positive-definite tangent stiffness exhibiting microdamage, show- 
ing the elastic ( W ) and irrecoverable ( S ) portions. 
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Figure 3: Schematic showing the transition form a continuum to a cohesive zone due to the initiation 

of macroscopic cracks. The essential, constitutive variables switch from stress and strain to traction and 
separation. 





(c) Mode II matrix fracture. 

Figure 4: Triangular traction versus separation which dictates the behavior of cohesive cracks embedded in 
the continuum. The total area under the traction-separation law represents the material fracture toughness 
G J mC . The area above the unloading line for a given traction-separation state is the strain energy release 
rate G 3 m . 
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Figure 5: Example problem used to demonstrate mesh objectivity of EST. One quarter of a 200 mm x 100 
mm containing a central hole with a radius of 5 mm is loaded in tension with symmetric boundary conditions 
on the bottom and left boundaries. 


L 


(a) 0.57V 


(b) 0.2rv 





(c) 0.1r h . 


(d) 0.0477, . 


Figure 6: Four mesh densities used to demonstrate the mesh objectivity of EST. 


21 of 40 


American Institute of Aeronautics and Astronautics 



Figure 7: Reaction stress normalized by critical axial strain times axial Young’s modulus verse applied 

displacement normalized by hole radius for four different mesh densities. 



(a) 0.5r/j,, a = 0.50. (b) 0.2r^,, a = (c) O.lr^,, a = 0.37. (d) 0.04r^,, a = 

0.42. 0.36. 


Figure 8: Contours of the reduced microdamage ISV 5 r , normalized by the maximum S r obtained from all 
simulations, immediately prior to failure initiation for four different mesh densities. 
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(a) 0.5r^, cr = 0.80. 



(c) O.lr^, cr = 0.82. 



(e) 0.5 rhi cr = 0.030. 



(g) O.lr/j,, a = 0.042. 



(b) 0.2r^, cr — 0.81. 



(d) 0.04r/j,, a = 0.82. 



(f) 0.2rh, a = 0.018. 



(h) 0.04r^,, a = 0.242. 


Figure 9: Contours of the transverse degradation parameter indicative of the transverse crack path in 
the specimens. The contours (a-d) are presented at the ultimate load achieved by the specimens and (e-h) 
after the specimens have lost their load carrying capability. 
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Figure 10: Geometry, boundary conditions, and strain gage (Sg) locations of CNPs tested at NASA LaRC. 39 
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Figure 11: FEM mesh used to simulate tensile loading of CNPs. 
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(a) Shear microdamage function obtained from (b) Transverse tension and compression micro- 
[±45° ]s angle-ply laminate. damage functions obtained by scaling data for 

AS4/3502 in Ref. 36. 

Figure 12: Microdamage functions for T800/3900-2 used in FEM models. 



Figure 13: Applied load versus displacement of a 4” section for Laminate- 1. 
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P (lbf.) P (lbf.) 




(a) Sg-1. (b) Sg-2. 




(c) Sg-3. (d) Sg-4. 

Figure 14: Applied load versus local strain for Laminate-1 


Notch 



Damage Path (Splitting) 


Figure 15: C-Scan of failed Laminate- 1 specimen. 39 
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(a) P = 8,210 lbf. 



Figure 16: 


(b) P = 16,400 lbf. 

Normalized matrix microdamage contour s ^ r ax 


in Laminate- 1. 
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(a) P = 8,210 lbf. 



(b) P = 16,400 lbf. 

Figure 17: Matrix shear failure degradation DJj in Laminate- 1 
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16000 



Figure 18: Applied load versus displacement of a 4” section for Laminate- 2. 
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P (lbf.) P (lbf.) 




(a) Sg-1. (b) Sg-2. 




(c) Sg-3. (d) Sg-4. 

Figure 19: Applied load versus local strain for Laminate-2. 



Debonding of 45° ply 


Figure 20: Photograph of failed Laminate-2 specimen. 39 
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(b) 0° Layer. 



(c) -45° Layer. (d) 90° Layer. 


Figure 21: Normalized matrix microdamage contour J^ ax in Laminate-2 just prior to first axial failure 

initiation P = 8,640 lbf. 
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(a) 45° Layer. (b) 0° Layer. 



(c) -45° Layer. (d) 90° Layer. 


Figure 22: Fiber failure degradation D{ in Laminate-2 at ultimate load P = 15,300 lbf (magnified view of 
region near notch). 


32 of 40 


American Institute of Aeronautics and Astronautics 





Field-2 


SNEG, (fraction = -1.0) 
(Avg : 75%) 


S 


+1.000e+00 
+9.167e-01 
+8.333e-01 
+7. 500e-01 
+6.667e-01 
+5.833e-01 
+5.000e-01 
+4. 167e-01 
+3.333e-01 
+2. 500e-01 
+1.667e-01 
+8.333e-02 
+0.000e+00 




(a) 45° Layer. 



(b) 0° Layer. 



(c) -45° Layer. (d) 90° Layer. 

Figure 23: Transverse matrix failure degradation D™ in Laminate-2 at ultimate load P = 15,300 lbf. 
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(a) 45° Layer. 



(b) 0° Layer. 



(c) -45° Layer. 

Figure 24: Shear matrix failure degradation D y\ 


(d) 90° Layer. 

in Laminate- 2 at ultimate load P = 15,300 lbf. 
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(a) 45° Layer. 



(b) 0° Layer. 



(c) -45° Layer. 



(d) 90° Layer. 


Figure 25: Normalized matrix microdamage contour J^ ax in Laminate-2 after specimen has lost load 

carrying capability. 
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(a) 45° Layer. 



(b) 0° Layer. 



(c) -45° Layer. 


(d) 90° Layer. 


Figure 26: Fiber failure degradation D{ in Laminate-2 after specimen has lost load carrying capability. 
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(a) 45° Layer. 


(b) 0° Layer. 



(c) -45° Layer. (d) 90° Layer. 


Figure 27: Transverse matrix failure degradation D ™ in Laminate-2 after specimen has lost load carrying 
capability. 
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(a) 45° Layer. 


(b) 0° Layer. 



(c) -45° Layer. (d) 90° Layer. 


Figure 28: Shear matrix failure degradation DJ} in Laminate-2 after specimen has lost load carrying 

capability. 
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